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Abstract. - The fluctuations of the work done by an external Gaussian random force on a 
harmonic oscillator that is also in contact with a thermal bath is studied. We have obtained 
the exact large deviation function as well as the complete asymptotic forms of the probability 
density function. The distribution of the work done are found to be non-Gaussian. The steady 
state fluctuation theorem holds only if the ratio of the variances, of the external random forcing 
and the thermal noise respectively, is less than 1/3. On the other hand, the transient fluctuation 
theorem holds (asymptotically) for all the values of that ratio. The theoretical asymptotic forms 
of the probability density function are in very good agreement with the numerics as well as with 
an experiment. 



One of the most fundamental and important problems 
in the nonequilibrium physics is to understand fluctua- 
tions. In this context, the so-called fluctuation theorem 
(FT) has generated a lot of interest. The FT was found 
first for the phase space contraction in dynamical sys- 
tems [l][2j and later for a certain "action functional" in 
stochastic systems [3[|4] — these quantities are generally 
referred to as the "entropy production". Subsequently, 
there has been an increased interest in the FTs for var- 
ious physical quantities such as work, power flux, heat 
flow, total entropy, etc. [5||7] — because, in the absence of 
a general framework for nonequilibrium phenomena, the 
FTs seem to be providing an unifying picture for a va- 
riety of nonequilibrium systems. The so-called Jarzynski 
equality [S], Crooks relation f9l, and Hatano-Sasa iden- 
tity [To] are closely related to the FT. In the linear re- 
sponse regime, the FT leads to the Green-Kubo formula 
and the Onsager reciprocity relations 4pT . However, the 
FT is more general, as it also describes fluctuations in the 
nonlinear regime arbitrarily far from the equilibrium. 

The FT relates the positive and the negative fluctua- 
tions of a certain time-integrated physical quantity Wr = 
W{t) dt, during a nonequilibrium process, according to: 
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P{Wr = -Wt) 



(1) 



where P{Wr — ±wt) is the probability density function 
(PDF) of the physical quantity WV to have a value ±wt. 



In fact, depending on the choice of the initial ensemble, 
there are two kinds of FTs: the transient fluctuation the- 
orem (TFT) — in which the system at t = is in equilib- 
rium, and the steady state fluctuation theorem (SSFT) — 
in which the quantity WV is computed in a time interval 
T in the nonequilibrium steady state. Usually, the TFT 
is stated for a finite r, i.e., without the limit t — >■ oo in 
eq. Naively, one would expect the TFT and the SSFT 
to become equivalent in the t — > oo limit. However, this 
is not always correct. 

There have been several experimental tests of the FT 
and related results, in diverse systems such as a colloidal 
particle in a changing optical trap 
tal electroconvection 



12 ■ 14 , liquid crys- 



15 , fluidized granular medium 16 



cellar gel 20 , harmonic oscillator 21 , self-propelled po- 



lar particle 22 , wave Turbulence 23 , and a gravitational 
wave detector [24'. A recent review of the experimental 



electrical circuits 17 , RNA stretching T8p9 , sheared mi 



applications of the FTs may be found in ref. 25 . Interpre- 



tation of experimental findings are not always easy as the 
FTs are governed by the atypical fluctuations that corre- 
spond to the tails of the probability distributions — and 
in an experiment in a finite time, it is often hard to ac- 
quire enough of the rare events to produce the tail of the 
distribution accurately. Therefore, it is very important to 
have exact theoretical predictions. 

Theoretical investigations of the work FTs so far have 
been mostly limited to the systems describe by linear 
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Langevin equations with a Gaussian white thermal noise 
and driven out of equihbrium by an external deterministic 
force. In such cases [6j , the distributions of the work done 
by the external force are Gaussian and hence the work FTs 
hold somewhat trivially. On the contrary, the distribu- 
tions of the work done by an external Gaussian stochastic 
force have been found to be non-Gaussian in recent exper- 
iments on systems coupled to a thermal bath and driven 
out of equilibrium by an external random force 26 . Mo- 
tivated by these experiments, in this Letter, we address 
the important question regarding the role of the external 
stochastic forcing on the work fluctuations. 

We consider one of the most basic physical systems, 
namely, the harmonic oscillator. We investigate the fluctu- 
ations of the work done by an externally applied Gaussian 
random force on a harmonic oscillator that is also in con- 
tact with a thermal bath. The displacement x{t) of the 
harmonic oscillator from its mean position is described by 
the Langevin equation 



dt^ 



+ l^ + kx^CTit) + Mt), (2) 

where m is the mass, 7 is the viscous drag coefficient and 
k is the spring constant. The interaction with the ther- 
mal bath is modeled by a Gaussian white noise Ct(^) with 
zero-mean (Ct(*)) = 0. The externally applied force /o(t) 
is again a Gaussian random variable with {fa(t)) — 0, 
and and /o arc uncorrelated. Equation ([2| is asym- 
metric in and /o — the fluctuation-dissipation theo- 
rem relates the thermal fluctuation to the viscous drag as 
(CT(s)CT(i)> ^ 2DS{s - t) where D = 7/cbT with T be- 
ing the temperature of the bath and fc^ being the Boltz- 
mann constant, whereas the fluctuation of the external 
force (/o(s)/o(i)) — [5 fo)'^5{s — t) is independent of 7. As 
it turns out, the only relevant parameter is 



{6h 



2D 



where {x ) and 



/cq 



— 1, and a e (0,00), (3) 



are the variance of x in the steady- 



state (for /o 7^ 0) and in equilibrium (for /q = 0) respec- 
tively. 

The quantity of interest is the work done by the external 
random force /o(t) on the harmonic oscillator in a time 
interval r, in the nonequilibrium steady state. This is 
given (in units of kgT) by 



Wr 



1 



n , -.dX , 
hit)-dt. 



(4) 



with the initial condition (at t — Q) drawn from the steady 
state distribution. Evidently, WV is a fluctuating quantity 
whose value depends on the initial condition, the trajecto- 
ries of thermal noise {C,T{t) '■ ^ <t < t} and the external 
random force {/o(i) : < t < t}, during any particular 
realization. 

It is clear from eq. ([2| that both the displacement x 
and the velocity v — dx/dt depend linearly on the ther- 
mal noise and the external random force. Therefore, the 



distribution of the phase space variables (x, v) is a Gaus- 
sian whose covariance matrix can be easily evaluated from 
eq. ([2]). However, due to the nonlinear dependence of the 
work given by eq. Q, on the thermal noise and the exter- 
nal random forcing, the PDF P(WV) is not expected to be 
Gaussian — although for any fixed realizations of {fo{t)} 
the work fluctuation would be Gaussian. Nonetheless, one 



expects the large deviation form 27 



P{Wr = wt/t^) - e^^/^^) for r > r. 
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(5) 



where Tj = m/"f is the viscous relaxation time and h{w) 
is the large deviation function (LDF), which is defined by 



h^w) = lim 



1 



lnP{Wr = wt/t^). (6) 



The FT as given by eq. ([I]) is equivalent to the symmetry 
relation 



h{w) — h(—w) — w. 



(7) 



Our aim is to obtain the LDF h{w) exactly, as well as the 
complete asymptotic form of the PDF P{Wr). 

We begin by considering the characteristic function 



-\W-r\ 



dWre-^^^P{Wr) = Z{X,t), (8) 



where (• • • ) denotes an average over the histories of the 
thermal noise and the random forcing as well as the 
initial condition. The restricted characteristic function 
Z(X,x,v,t\xo,vo) — where the expectation is taken over 
all trajectories of the system that evolve from a given ini- 
tial configuration {xq,vq) to a given final configuration 
{x,v) in time r — satisfies the Fokkcr-Planck equation 
— C\\ Z{X, x, V, t\xq, vq) = with the initial condition 
Z(A, X, w, 0\xo, Vq) = S{x — xo)S{v — Vq), whcrc the Fokker- 
Planck operator is given by 



—x + —{l + 2aX)v 
m m 

aX'^j^ 2 , 7 



dv 



dx 



D 



v' + ^{l + aX). (9) 



The solution of the Fokker-Planck equation can be for- 
mally expressed in the eigenbases of the operator C\ and 
the large r behavior is dominated by the term having the 
largest eigenvalue. Thus, for large t, 

Z{X,x,v,t\xq,vq) - x(a;o,wo,A)*(a;,?;,A)e^^(^), (10) 

where ^'(a;, w, A) is the eigenfunction corresponding to the 
largest eigenvalue ^{X) and x(a;o, vq, X) is the projection of 
the initial state onto the eigenstate corresponding to the 
eigenvalue ^(A). To calculate these functions, we follow an 
approach that was used recently to compute the fluctua- 
tions of the heat transport across a harmonic chain l28| . 
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Skipping details 1291, we find that 



^,{\) = ^[l-r,{X)], ,7(A) = v/l + 4aA(l-A), (11) 



*(x,i>,A) 



2TT{l + a)D 
and x{xo,vo,X) = exp[-B^{X)E{xo,vo)], 



exp[-B+{X)E{x,v)], (12) 
(13) 



where 



and 



B±iX) 



7[77(A)±(1 + 2qA)] 
2(l + a)D 



E{x, v) = -^kx^ 



(14) 



(15) 



is the total energy of the harmonic oscillator. Note from 



eq. (11 1 that the largest eigenvalue satisfies the symmetry 
relation /i(A) — /i(l — A), even though L\ and its adjoint 



do not possess the symmetry = Ci-\. 
Using the explicit forms of eqs. Q and ([ll])~( 13 ), the 
eigenvalue equation C\^{x, v, A) = /i(A)5'(a;, v, A) and the 
normalization x(a;, u, A)^'(a;, u, A) da; = 1 can 

be indeed verified. Moreover, ^(0) = and x{xo,vo,0) = 
1, which is expected — since eq. ^ for A = 0, corresponds 
to the Fokker-Planck operator of the phase space variables, 
and hence the steady state distribution Z{X = 0, x, u, r — > 
oo\xo,Vq) must be independent of the initial condition and 
T. The steady state distribution of (x, v) is given by Z(X = 
0,x,v,T cx3|xo, fo) — ^'(a;, w, 0). 



Now, substituting Eqs. (12) and (131 in eq. (10), then 



averaging over the initial variables {xo,vq) with respect 
to ^(0^0,^0,0) and integrating over the final variables 
{x,v), we find the characteristic function that is defined 
by eq. as 

Z(A,T)^g(A)e-^(^), (16) 



where /^(A) is given by eq. ( 11 ) and 



3(A) 



2r,(A) 



1 + 77(A) -2q:A l + r7(A)+2aA 



(17) 



The first factor in the above equation is due to the aver- 
aging over the initial conditions with respect to the the 
steady state distribution and the second factor is due to 
the integrating out of the final degrees of freedom. 

The PDF of the work done is related to its characteristic 
function by the inverse Fourier transform 



1 

27ri 



-\-ioo 



Z{X,t) e 



dX, 



(18) 



where the integration is done along the imaginary axis 
(vertical contour through the origin) in the complex A 
plane. The large r (^ t^) behavior of P{Wr) can be ob- 
tained from the saddle point approximation of the above 
integral while using the asymptotic form of Z(X, t) given 



by eq. (16). We note that 7?(A), given in eq. (11), has two 
branch points on the real A line at 



l± W1 + - 
a 



(19) 



■\/4q;(A+ — A)(A — A_) . Outside the interval 



on the real A line, 7y(A) is imaginary. However, 



as 77(A) 

[A-,A4 

Z{X, t) must be real for real values of A, if the integral 
in eq. (|8| converges. Therefore, analytical continuation of 
Z{X) to the real A is allowed only within the range A_ < 
A < A+ — for which [77(A)]^ > 0, and hence, /i(A) is real 
and analytic. In fact, in the whole complex A plane, 77(A) 
is real only for A in the real interval [A_,A+]. Therefore, 
we expect the saddle to be also in that interval. 

Now, in the expression of g{X) given by eq. ( |17p , the de- 
nominator of the second factor is positive for A € (A_, A+) 
for all a G (0,oo). Hence, the second factor of g{X) is an- 
alytic in the interval (A_,A+). On the other hand, the 



analytic properties of the first factor in eq. (17), depends 
on the value of the parameter a. 

As long as a < 1/3, the denominator of the first factor 
is positive for A e (A_,A+). Therefore, in this case g{X) 
is analytic in (A_, A+) and hence can be neglected in the 
saddle-point calculation as a subleading contribution. The 
saddle-point calculation with Z{X, r) ~ e^^^^'^ relates /i(A) 
to the LDF h{w) of eq. ([5|, by the Legendre transform 



h[w) = T^^[X*) + A*w, 



(20) 



In this case, the symmetry relation of the LDF as given 
by eq. ([t]), follows directly from the symmetry /i(A) = 



/i(l — A). The solution of the condition - 
gives the saddle point A* in terms of w as 



yM'(A*) 



X*{w) 



(21) 



We now consider the case a > 1/3. In this case, due to 



the first factor in eq. (17), g{X) possesses a pole at 



Ao = 



1 



(22) 



and A_ < < Ao < A+. Now, g(X) is negative for A > 
Aq. However, (?(A) must be non- negative for any real A, if 
the integral in eq. (Isl) exists. Therefore, now the allowed 



range of real A shrinks to (A_, Aq). It follows from eq. (21 ) 
that X*(yj) is a monotonically decreasing function of w, 
and X*(w =foo) A±. Note that A* e (A_,A+) as 
expected. For any given a, as w decreases from -fcx) to 
— cx), the saddle point X^iw) moves unidirectionally from 
A_ to A+. Thus, for sufficiently large w, we have A_ < 
A* < Aq. In such situation, the contour of integration can 
be deformed smoothly through the saddle point A*, and 
therefore, the LDF is still given by h(w) — T^yuiX*) + A*w. 
However, as one decreases w, at some particular value w — 
w* , the saddle-point hits the singularity. For w < w* , 
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Fig. 1: (Color online). Plot of the large deviation function 
h{w) for a = 1/5 (red), 1/3 (green), 1 (blue), 2 (brown), and 
3 (magenta). 



we then have < Aq < A*. In this case, the leading 
contribution comes essentially from the pole [29], which 
yields h^w) = r^//(Ao) + XqW. Using X*{w*) = Aq and 
/i'(A*) + w = 0, it is easy to check that h{w) and its 
derivative h'{w) are continuous aX w — w* . For a = 1/3, 
we have Aq = A+ = 3/2. Since, A* — > A+, only when w — > 
— oo, for any finite w we again have h{w) — T-y^(A*) + A*w, 
i.e., w* — —oo. 

Let us express the LDF h{w), defined by eq. (|6]), explic- 
itly in terms of w and a. We find that, for a < 1/3: 




h{w) — hi{w), 



hi{w) for w > w* 



h2[w) for w <w* ^ 
where hi{w) and h2{'w) are given by 



(23) 



(24) 



hi{w) 
h2{w) 



— \fvp- 
2w 



1 

1 + - 
a 



1+a 1+a 
and w* is found by solving X*{w*) 



Aq, as 



w 



a{a — 3) 
3a- 1 ■ 



(25) 
(26) 

(27) 



Figure [T] displays the LDF for various a. 

From the above expressions, it is now straightforward 
to check the validity of the work SSFT. For a < 1/3, we 
get h{w) — h{—w) = w, which implies that the SSFT is 
satisfied. On the other hand, for 1/3 < a < 3, we get 
h{w) — h{—w) = w only for w* < w < —w*. For a > 3, 
the symmetry relation ([t]) is not satisfied for any w. For 




Fig. 2: (Color online). Asymmetry function for various a. 



example, for a = 3, we get w* = and h{w) — h(—w) = 
1 + ui — y^l + {w^/3) . Figure [i] displays the asymmetry 
function p{w) = h(w) — h{—w) for several values of a. 
We can also obtain the the complete asymptotic form 



of the PDF of the work fluctuations. Skipping details 29 
we find that for a < 1/3: 



K{w) 



^{T/T-,)hi{w) 



see eq. ( 29 I 



P{Wr = wr/r-y) 
and for a > 1 /3 
where 

K{w) = a3/2(l + l/«)3/4(u,2 + a)-5/4 

r 1 ^1 

X 1 + {w + a)X*{w) ~ hi{w) 

X \{h2{w) - hi{w)} + {w- a){X*{w) - Aq} 



(28) 



and 



= hn. [(A - Ao)5(A)] = -^^T^ 
A->-Ao 8q;^(1 + a) 



(30) 



(31) 



We compare the above asymptotic forms of the PDF 
with the results obtained from the numerical simulations 
of the Langevin equation As seen from fig. [sj the 



agreements are extremely good, even for r/r-y 



10. 



We also compare our analytical results against an exper- 
iment that was carried out on an atomic-force microscopy 
cantilever and reported in ref. 26 . The dynamics of the 



micro-cantilever tip is described by eq. ([2| with the viscous 
relaxation time r-y — 632 fis. We find very good agreement 
between the theory and the experiment, which is shown 
infig.|4] 
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P(Wr = wt/t^) 



^{T/T^)hi{w) 



K[w) 



sgn(t«* 



w) g^i 



^Jh2{w) ~ hi{w) 



sgn(it;* — w) 



eTk(^y/T[h2{w) - hi{w)]j - e{w* - w) 



(29) 
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Fig. 3: (Color online). P{Wt) against the scaled variable w — {t /T.y)~^Wr for (t/t^) = 10 and (a) a = 1/3, (b) a = 3 and 
(c) a = 20 respectively. The points (blue) are obtained from numerical simulation, and the thick dashed lines (red) plot the 
analytical asymptotic forms given by eq. (281 for (a) and eq. (291 for (b) and (c). The vertical thin dashed lines in (b) and (c) 
mark the position of w — w* . 
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Fig. 4: (Color online). P{Wt) against the scaled variable w — (j lT^')~^Wr for {t/t^) = 7.88..., and (a) a = 3.03 and (b) 
a = 18.66 respectively. The points (blue) are from the micro-cantilever experiment reported in ref. [26], and the thick dashed 
lines (red) plot the analytical asymptotic form given by eq. (291. The vertical thin dashed lines mark the position oi w = w* . 
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So far, we have considered the fluctuations of Wr in the 
nonequilibrium steady state, as we have averaged over the 
initial conditions in eq. ( |10[ ) with respect to the steady 
state distribution ^'(xo,wo,0) to arrive at eq. ( [l6| . Let 
us now examine how the nature of the initial state af- 
fects the results. We recall that the singular part of 5(A), 
i.e., the first factor in eq. (17) comes from the averag- 
ing of eq. (10) with respect to the steady state distri- 



bution of the initial state. Without the averaging, for 
any given initial configuration (xo,vo), the resulting pref- 
actor of e'^''^'^^ remains analytic throughout the interval 
(A_ , A+) , and hence can be neglected from the saddle point 
calculation as the subleading contribution. Therefore, the 
FT for a fixed initial condition is always satisfied, as the 
LDF is given by eq. (23) for all a S (0,oo). If the initial 
state at r = is chosen from equilibrium — i.e., the av- 
erage in eq. ( 10 ) is taken with respect to the Boltzmann 
weight cx exp[—E{xo,vo)/{kBT)] — then the first factor 
in eq. ^ is replaced by 2(1 + a)/[l + r?(A) + 2a(l - A)]. 



In that case even g{X) satisfies the symmetry relation 
g{X) = g{l — A). It is easy to see that, now g{X) re- 
mains analytic in (A_, A+) for any a. Therefore, the LDF 
in this case is again given by eq. (23) for all a € (0, 00). 
Consequently, the TFT is satisfied (as t — > 00) for all 
a S (0, 00). 

In conclusion, we have studied the work fluctuations of a 
harmonic oscillator coupled to a thermal bath and driven 
out of equilibrium by an external Gaussian random force. 
We have found that the SSFT holds only for weak forc- 
ing, whereas the TFT (with t — >■ 00) holds for all forcing. 
More importantly, we have analytically obtained the exact 
LDF as well as the complete asymptotic forms of the PDF 
of the work fluctuations, and quite interestingly, they are 
independent of the spring constant of the harmonic oscil- 
lator. However, while the LDFs are same for both fc ^ 
and k = cases, the complete asymptotic form of the 
PDFs are different in the two cases. Therefore, A: — 
limit of PDF (which is anyway independent of k) is not 
same as the PDF in the A: = case. The nature of the 
work fluctuation is found to be non-Gaussian. These ex- 
act results should have broad and important applications, 
as the harmonic oscillator is ubiquitous in nature. For ex- 
ample, many nanomechanical and biological systems are 
essentially described by a harmonic oscillator and the re- 
sults of this Letter are expected to be useful there. 
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